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ABSTRACT 

A time-series technique is presented for identifying the 
dynamic characteristics of the human operator in manual control 
tasks from relatively short records of experimental data. Con- 
trol of system excitation signals used in the identification is 
not required. The approach is a multi-channel identification 
technique for modeling multi-input /multi-output situations . 
The method presented includes statistical tests for validity, 
is designed for digital computation, and yields estimates for 
the frequency responses of the human operator. A comprehensive 
relative power analysis may also be performed for validated 
models. This method is applied to several sets of experimental 
data; the results are discussed and shown to compare favorably 
with previous research findings. New results are also 
presented for a multi-input task that has not been previously 
modeled to demonstrate the strengths of the method. 
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J_. INTRODUCTION 

A pilot model is a mathematical expression which balances simplicity 
of mathematical structure with observed empirical reality according to the 
purpose for which it is used. A key question always facing the aviation 
community has been how to develop and use these models in order to specify, 
design, and evaluate piloted systems 1 so that they provide efficient, pro.- 
ven ^performance while admitting the pilot "symbiotically" into the control 
loop . The successes of describing function and optimal control models in 
meeting this objective are well known 3 , but the identification of these 
models is hindered by an dependence on long data records, a priori parame- 
ter knowledge, and a precisely controlled experimental environment. 

A time series approach to pilot modeling, introduced ten years ago^, 
initially appeared as just another "technique"; but recent applications of 
time series analysis to complex multi-channel tasks 3 indicate that this 
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approach may work well on relatively short data records with little or no a 
priori parameter knowledge. Moreover, the process of modeling provides a 
unifying mathematical "framework" relating recent research in closed-loop 
multi-channel identification theory to actual laboratory or flight test 
data records of relatively short duration. The "framework" includes estab- 
lishing model existence, applying a proven identification technique, vali- 
dating the resulting model, and analyzing model properties relative to 
model purpose. 

Early researchers using time series to model manual control behavior 
recognized that obtaining single or multi-channel pilot models is a doubly 
formidable task because of the adaptive nature of the pilot and because of 
the inherent loop closures in the overall system 6 . Shinners and Agarwal , 
in their pioneering work for single-input, single-output (dual-channel) sys- 
tems, found that simple discrete transfer functions adequately described 
pilot manual control output in compensatory and pursuit tasks but did not 
consider the theoretical question of model existence or stability. The 
work of Goto, based on the theoretical methods of Akaike and Whittle , 
considered model "existence" questions for a two subsystem closed-loop 
structure 10 , but these methods assume that the autocorrelation statistics 
for the process are known a priori. 

The purpose of this paper is to provide a unifying framework for time 
series modeling by deriving the specific theoretical and experimental con- 
ditions required for model existence and uniqueness, to apply an identifi- 
cation algorithm which guarantees stability and does not require a priori 
statistical information, and to demonstrate the application of this iden- 
tification process in case studies. The derivation of existence conditions 
is applicable to a three subsystem closed-loop structure which contains the 
two subsystem results of Goto as a special case. The derived identifica- 
tion algorithm is called "Normalized Predictive Deconvolution", NPD, and is 
a generalization of the Levinson-Wiggins-Robinson algorithm and the 
multi- channel Maximum Entropy Spectral Estimation algorithm . 


2 . THE MODEL 

The pilot-as-controller discrete linear model is shown as part of a 
three subsystem structure in Figure 1. The double lines represent vector 
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precesses from three subsystems: the vehicle, the pilot, and the flight 

control system. Autoregressive (Markov) noise is added to each subsystem 

13 

to represent a physical disturbance ; that is, injected noise is a linear 
sum of past values plus an i.i.d. discrete "shock" or "pulse". Mathemati- 
cally this representation may be concisely represented by 


X(t) = G(z)X(t) + Y(t) 

where X(t) is a joint process vector partitioned into subsystems as 

X(t) = |f T (t),6 T (t),e T (t)| 


( 1 ) 

( 2 ) 


G(z) is a matrix of transfer functions in terms of the shift operator "z" 
which may also be partitioned into a general form given by 


G(z) 


0 G 12 (z) 

G 2i<z)= G f (z) 0 

G 31 (z> G 32 (z)==G a (z) 


G l 3 (z)= G p ( z) 

g 23 (z) 

0 


(3) 


The injected noise, V (t), is assumed both autoregressive of finite order 
"L" and uncorrelated between subsystems. Thus, it may be represented by 
the block diagonal form 

f(t) = C(L,z)¥(t) + p(t) (4) 


C(L,z) 


2 C 11 k z 
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-k 
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!, C 33,k z 


-k 


(5) 


¥(t) = R T (t) , V T (t) , W T (t) 


( 6 ) 


4 



p(t) 


r T (t) , v T (t) , w T (t) 


i.i.d. 


(7) 


The individual elements in Equation (3), in contrast to the finite 
order assumption for the noise representation, are expressible either as a 
ratio of discrete polynomials (transfer function) or as an infinite 
sequence in the delay operator z (pulse response). Thus, between subsys- 
tems "i" and "j". 


G ij (z) = ^ij.k z 


-k 


( 8 ) 


If the infinite sequence of Equation (8) is truncated at order "M", an 
approximation to the mathematical system of Equation (1) results which will 
be referred to as the joint autoregressive representation (JAR). The trun- 
cated elements of G(z) are given by 


V M > 2) " J/W 2 


-k 


JAR 


(9) 


By combining Equations (1), (4), and (9) the JAR may be written as 

X(t) = G(M,z) X(t) + «(z) p(t) (10) 


« _1 (z) = I - C(L,z) 


(ID 


C il (L,z) = Z C 


-k 




( 12 ) 
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The joint innovations representation , JIR, is obtained by multiply- 
ing Equation (10) by Equation (11) and solving for X(t): 


X(t) *= A(M,z) X(t) + p(t) 


M 


-k 


A(M,z) = I A, k z s 
k=l 


C(L,z) + G(M,z) - C(L,z)G(M,z) 


(13) 

(14) 


The block diagonal form of Equations (2) and (5) .is now taken into 
account in the relationship between the JAR and JIR. Denoting each 
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subsystem of X(t) by subscript "i", Equations (13) and (14) are equivalent 
to 


X (t) - C (L,z) X, (t) + L 

J-l 


I - C 1± (L t z) 


G lj (M,z) X^t) + Pl (t) (15) 


By comparing Equations (14) and (15) one obtains 


Cii( L ,z) = A ljL (L,z) 


(16) 


“ C 14 (L,z) G, .(M.z) + A., .,(M,z) ; i * j 


ii 


ij 


ij' 


(17) 


The JIR described by Equation (13) may also be put into the form 


X(t) = T(M,z) p(t) 


r(M,z) = 


M 


I - 


^ k Z 
c=l 


-k 


-1 


(18) 

(19) 


The autocovariance matrix is found by post multiplying Equation (18) 


by the transpose of X(t) and taking the expected value: 


R (0) = E 
xx 


X(t)X T (t) = r(M,z) P(0) r*(M,z) 


where 


P(0) - E 


p(t) p (t) 


( 20 ) 


( 21 ) 


The power spectral density of this process'* is 
* xx (io) = |r(M,z) A P(0)r*(M,z) 


icoA 

z = e J 


which has the property 


$ (z) = $ r *'(z^)=$ (z) 

xx xx xx 


( 22 ) 


(23) 


An approximation to the frequency response between variables "i” and "j" 
may be found using 


G ± j(w) ~ |G iJ (M,z)| 




(24) 


z ■ e 
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If P(0) is diagonal, the relative power in state "i" is defined as 


n 


P,,(£o) = E r .(u>) 
j=l J 




( 25 ) 


and the noise power contribution to channel "i" from the noise source in 
channel "j" is 


q lJ ( “ ) ‘ 4 V 0>r U < .“ ) ^ 


(26) 


Thus it is shown how the JIR representation of Equation (13) may be 
transformed into the JAR representation of Equations (9) through (12) using 
the recursions of Equations (16) and (17). Once validated, the properties 
of the identified model may be analyzed using Equations (20) through (26). 
There must be assurance, however, that these model representations exist in 
theory, and this topic is addressed in the next section. 


_3 . THE EXISTENCE QUESTION 

The primary factors in the determination of an acceptable pilot model 
are suitable experimental conditions, the assumed model structure, and the 
identification technique. Since the harm done by a faulty experiment, simu- 
lation, or flight test permanently voids the data, the conditions required 
for a unique and valid model are very important. 

THEOREM 1: The JIR of Eqn.(18) is unique, and there is a unique map- 
ping between the JAR of Eqn. (10) and the JIR of Eqn.(13) providing 
Eqn.(23) holds for the spectral density and providing there is a delay in 
every path of Figure 1. 

For proof see the Appendix. 

THEOREM 2: Given that the transfer matrix f(z) has been identified 
from realization set {X(t)| t<N } generated by T(z), necessary conditions 
for 


lim f(z)=r(z) (27) 

M,N-m» 

are 


(1) The joint process X(t) is full rank 
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(2) There is a unique factorization 


$ xx (u) “ r(z)UA 1/2 (UA 1/2 ) T r*(z) (28) 

z=e^ w 

A > 0 and U Unitary (29) 

For proof see the Appendix. 

I 

The practical implications of these theorems for flight simulations 
and flight tests are that sufficient noise sources be used to excite the 
vector process X(t), that there should be no feedforward paths which 
violate the requirement for a delay in each loop, and that no anticipatory 
loops are closed by the pilot for the same reason. Although some identifi- 
cation schemes allow correlated noise inputs 15 , there is no way to distin- 
guish them from feedforwards and/or anticipation. If validation tests, 
however, indicate a positive definite and diagonal autocorrelation matrix 
for the noise inputs, then there is evidence that a sufficient condition 
has been met for uniqueness. 

To summarize, the design or test engineer should assure 

(1) sufficient noise excitation in measured channels; 

(2) pilot anticipation negligible (implies random or random appearing 
inputs; 

(3) physical delays exist in each channel, including feedforward, which are 
significant relative to sample time; 

(4) data realizations are not predominantly unstable or nonstationary; 

(5) validation checks include a whiteness test for the estimated noise 
realizations. 

4. MODEL IDENTIFICATION AND VALIDATION 

Given the conditions are met for model existence, an identification 
scheme is desired which identifies the JIR of Equation (13) from data real- 
ization set {X(t)| t<N }. It is especially important that the scheme be 
stable (identified parameters are bounded) and not be dependent on a priori 
knowledge of autocorrelation statistics. The identification technique 
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presented here is called Normalized Predictive Deconvolution (NPD), which 
acts directly on the data sets and results in a stable and parsimonious 
JIR. 

The basic principle of the NPD scheme follows that established by Wig- 
gins and Robinson^ who generalized Burg's^ recursion for single-channel 
systems by hypothesizing a set of backward predictors given by 

X(t) = B(M,z) X(t)+p'(t) 

M k 

B(M,z) = ZB ■ z K 

k-1 M * k 


(30) 

(31) 


p'(t) = 


T 

r'tt) , 


v' T (t) , 


T 

w' r (t) 


i.i.d. 


(32) 


T 

By post multiplying Equation (13) by X f (t-k) and Equation (30) by 
T 

X (t+k) , taking expected value, and expressing the result in a block 

matrix form, the "normal equations" of Reference (17) results 


-A < • « • A i A 
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i-B -B . 
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-B 


m, 1 
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Q w (m) 
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0 

V m) 


(33) 


T(m) 


R (0) 
xx 


R (~m) 
xx 


R (m) 
xx 


R (0) 

XX 


(34) 


m 


Q w (m) = R„„(0)- Z A . R (~k) 


xx 


w ** 1 “ > 

_j m,k xx 


(35) 


Q B (m > " \x 


(°)- ZB R(k) 
k=1 m »k xx 


(36) 


In the NPD scheme the solution to the "normal equations" is recur- 
sively generated as order "m" is incremented without knowing the autocorre- 
lation matrices a priori. The top and bottom rows of Equation (33) are each 
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weighted with invertible forward and backward prediction scaling matrices 
S^(m) and Sg(m) so that 

K nrfl,i = S I 1(m+1 > Vl.i J 0 < i < ffl+1 (37) 

" X m+1,0 = S A^ (38) 

Vl.i = S B 1(m+1) Vl.i 5 0 < 1 (39) 


“Vl,0 " S B (40) 

To derive the forward recursion formula (the backward recursion fol- 
lows analagously) , the scaled bottom row of the "normal equations" is mul- 
tiplied by an arbitrary but invertible matrix and added to the top row of 
Equation (33). Next, the order is incremented from "m" to "m+1" and the 
scaled results are expressed in the form 

^m+1,1 *** ^m+l.m+l 

B m+1 ,m+l B m+1 ,m * 


T(m) = 


q F (m+i) 


0 ... Qg(m+1) 


(41) 


By matching the terms of Equation (41) with the previously obtained linear 
combination of rows the following recursion results: 


‘ S A 1(m+1) S A (m) 


V i"Y ><!„* w>. MH 


-1 


“I, 


(42) 


Q F (m+l) = S A 1 (m+l) S A (m) 


- 1 , 


Qp(m) - S A (m)e F (m+l)Q B (m)e B (m+l) 


(43) 


B 


mhl ,i 


Sg^nH-l) S B (m) 




(44) 
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Q B (m+l> = s" 1 S B (m) 


Q B (m)-S B 1 (m)e B (nH-l)Q F 1 (m)e F (nrfl) 


where 0 < i < m+1 in the above expressions, and where the forward and 
ward prediction error matrices are given by 


By defining 


m 


e (m+1) = R (m+1)- E A R (k) 

* xx k=i m,m+l-k xx 


m 


e (m+1) = R(-m-l)- E B , R (-k) 

B xx . j m,m+l-k xx 


p (m+1 ) - S“ 1 (m)e F (m+l)S B i (m) 
PaOb+I) = I “ p(m+l)p' I '(m+l) 
P R (nr+l) = I - p^(m+l)p(m+l) 


rT, 


B 


S A^ m+1 ^ = ^A^ m ^~ e F^ nl+ * (m)E B (m+l ) 
S (nr+1) = S R (m)-e R (m+l)S. 1 (m)e„(m+l) 


B 


gvux/ g v / v / F ' 


it may be shown using matrix algebra that 

i-l 


PA /2 (m+l) 


= s” 1 Cart-1) S A (m) 


P B /2 (nt+l) 


-1 


= Sg^m+l) S B (m) 
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m+1 ,i 


p 1 /: W> 

a 


-1 
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m+1 ,i 


P B /2 (m+l) 


-1 


B 


W, i- pT(n+1 1 <-) V-’VwW-l 


(45) 


back- 


(46) 


(47) 


(48) 

(49) 

(50) 

(51) 

(52) 


(53) 


(54) 


(55) 


(56) 
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If the scaling matrices of Equations (37) and (39) are chosen to be 
the "identity" matrix, then the classical Levinson-Wiggins-Robinson (LWR) 
algorithm of Reference (11) results in a normalized form. If the scaling 
matrices are chosen so that 




<4 /2 « 



(57) 


s ; l <-) 




-1 


(58) 


then Equation (48) defines the Partial Autocorrelation Coefficient (PAC) 12 
matrix. In addition, if the following approximations are used: 


where 


p(m+l) 



(m) 


-1 


Rg /Z (m) 



(59) 


V m > “ 


n m 

Z I„(m,t) I (m,t) 
t=m+l * F 


(60) 


N 


Rp B (m) = Z Ip(m,t) T (m,t-l) 
t=m+l 


(61) 


N 

®j(m) = Z Ig(m,t~l) Tg(m,t~l) 
t=m+l 


(62) 


m 


lp(m,t)=S A (m)I_(m,t) = X(t)- Z A X(t-k) 

at k = 2 


(63) 


i B (m,t)=S B (m)I B (m,t) = X(t)- Z B m m _ k X(t~k) 

k™ 1 


(64) 


then the multi-channel Maximum Entropy Spectral Estimation algorithm of 
Reference (12) is obtained. 
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Morf, Vieira, and Kailath have shown that there is a one-to-one 
correspondence between the PAC matrices defined above and the autocorrela- 
tion matrices for a joint stationary process; moreover, they show that the 
characterization theorem of stochastic processes assures PAC matrices with 
singular values less that unity. 

To determine the order "M" at which the above recursion is stopped, a 
variation of the multi-channel Akaike rule^, as modified by the recommen— 
dations of Kashyap , is presented here as the PAC selection criterion. 
This criterion assumes that, as the estimates for the PAC matrix elements 
becomes smaller, they become more random, thus causing the determinant to 
also become random. To balance this effect with a term sensitive to both 
order "m" and number of channels "n" , the following expression was chosen 
as the PAC selection rule: 

A 

Jp(®) “ N log |det p(m)| + m(n) log N (65) 


The order resulting in the "first" minimum value as order "m" is incre- 
mented is chosen for the JIR. 

Validation is accomplished by testing the forward innovations for 
whiteness. These residuals are estimated using Equation (62) and the 
matrix set 


E 


ij?(M, t )ip (M, t— k) 


t < N ; 0 < k 


( 66 ) 


which is visually tested for whiteness over a reasonable number of lags 
"k". Plots of JIR statistics vs actual statistics (if available) and time 
histories of actual vs predicted JIR data may also be used. 

Summarizing, a technique called Normalized Predictive Deconvolution 
has been presented to identify a stable JIR of Equation (13) without a 
priori knowledge of the process autocorrelation matrices shown in Equation 
(34). The algorithm is initialized at m=0 with 


V 0) " V 0) " R xx (0) 

where Equations (60) and (63) are used to approximate R xx (0)» The scaling 
ma ^ r ^ ces are then chosen, as in Equations (57) and (58) for example, then 
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by definition 



S A<°> 


(67) 


pJ /2 (o) 


-1 


- v o) 


( 68 ) 


and Equations (38) and (40) are used to find S and B _. 

m,0 m,0 

The PAC matrix p(l) is then computed from Equations (59) through (64), 
from which Pa( 1) anc * are ^ oun< ^ us i n 8 Equations (49) and (50). The 
new forward and backward predictors are determined from Equations (55) and 
(56) for m=l and finally the value of the PAC selection rule using Equation 
(65) is found. If desired the order is incremented and the process 
repeated. 

Once the JIR is identified the JAR may be determined using Equations 
(16) and (17). The model characteristics are then calculated using Equa- 
tions (22) through (26). Case studies which demonstrate the application of 
this identification process and analysis are presented next. 


5 . MODEL ANALYSIS ; CASE STUDIES 

In order to demonstrate the application of the JIR identification pro- 
cess on actual data sets a multi-channel "piloted" simulation was accom- 
plished in the Flight Simulation Laboratory at Purdue University. Three 
pilots performed lateral bank angle tracking tasks using aileron deflection 
inputs with and without rudder deflection inputs for assistance. 

In addition to obtaining the data sets, the goal of the simulation was 
to obtain subjective pilot ratings and comments for three vehicle confi- 
gurations. The configurations were representative of large aircraft with 

the dutch roll modes selected to yield level 1, 2, or 3 handling qualities 

21 

as currently in military specifications . Table 1 summarizes the dutch 
roll characteristics and the corresponding pilot ratings and comments 
obtained during the simulation. Approximately 25 seconds (500 points at a 
20 Hz sample rate) were used for modeling from each data run which was typ- 
ically 60 seconds long. 
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The pursuit display shown to the pilot for the three-channel simula- 
tion (channels were aileron error, aileron deflection, and rudder deflec- 
tion) is shown in Figure 2. For the two-channel simulation the "ball in 
the window" portion of the display was masked and no rudder inputs were 
allowed. Note from the ratings and comments in Table 1 that there is a 
considerable degradation for each configuration between the two-channel and 
the three-channel cases. This degradation is most severe for the level 3 
configuration where a lateral pilot induced oscillation (PIO) resulted when 
the pilots were allowed to use rudder inputs. 

The commanded bank angle disturbance was a second order autoregressive 
process given by 

W(t) = 1.975 W(t-l) - 0.977 W(t-2) + .003 w(t) (69) 
w(t) = i.i.d. normal (0,1) (70) 

The parameters of this process were experimentally determined before taking 
tracking data to provide a realistic and unpredictable tracking signal to 
the pilots. 

The JIR pilot model was identified using the NPD algorithm set up to 

provide the special case of the multi-channel Maximum Entropy Spectral 

12 

Estimation algorithm . The PAC order selection rule of Equation (65) con- 
sistently resulted in M=4 in Equation (13) except for the three-channel 
Configuration 3 where the order was M-7. Figure 3 illustrates the behavior 
of the PAC selection rule versus order for this case. 

A typical experimental versus identif ied-model time history for the 
rudder deflection signal is shown in Figure 4 for models identified from 
100, 200, and 500 points. The 100 point model used every fourth point of 
the data set between points 1 and 400; the 200 point model used every other 
point between points 1 and 400. Thus the final five seconds of the time 
history shows actual and predicted time histories which are independent of 
the modeling process. The 500 point model shows the best visual agreement 
between actual and predicted time histories. 1 

The top row of the "normal equations" from Equation (33) may be used 
to define the predicted autocorrelation matrix as a function of lag for the 
identified JIR. With aileron deflection and aileron error as channels 1 and 
2, respectively, the actual versus predicted autocorrelation matrix is 
shown in Figure 5 for the two-channel Configuration 3 , where the actual 
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value was estimated from the data sets using 


R 

xx 


(k) = E 


X(t)X T (t-k) 


( 71 ) 


The normalized residual matrix from Equations (63) and (66) is shown 
in Figure 6. Normalization implies that each element is divided by the 
square root of the products of the respective diagonal element magnitudes, 
or 


NORMALIZED (i, j) 


ELEMENT(i, 1) 


\ 


ELEMENT (i,i) 


ELEMENT (j, j) 


(72) 


The prediction capability demonstrated in Figures 4 through 6 was typical 
for all identified models and was used as a validation check for all confi- 
gurations. From these results it was assumed that the models passed the 
validation checks using experimental data. 

If a model passes a validation check, the relative power analysis 
described by Goto 5 may be accomplished. The total power (variance) in the 
pilot's aileron deflection signal, computed from Equation (25), versus fre- 
quency for each two-channel configuration may be seen in Figure 7. Note 
that the power spectral density peak magnitude, in general, increases for 
configurations with higher (worse) pilot rating. Thus there is an indica- 
tion that pilot workload (as evidenced by power spectral density) increases 
across a portion of pilot bandwidth as pilot rating increases for different 

configurations. This is consistent with workload being correlated with 
22 

deflection rate . 

Using Equation (26) it is possible to calculate the amount of power 
due to the noise source in each channel. The noise contribution versus 
frequency for the aileron channel is shown in Figure 8 for the two-channel 
Configuration 3 (the other configurations showed similar results). Note 
that the command disturbance noise is the primary contributor to pilot 
aileron deflection at low frequencies (below 3 rad/sec) and pilot injected 
noise (remnant) is the primary contributor to pilot aileron deflection at 
the higher frequencies (above 6 rad/sec). The two-channel results are sum- 
marized In Table 2 and the three-channel results are summarized in Table 3. 
As expected, the error variance, or element (2,2) in columns 2 and 5 of 
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Table 3, increases both with pilot rating and with the added workload of 
the three-channel task (as measured by the spectral density). 

For the three-channel case studies, the total power in the pilot's 
aileron deflection signal for each configuration is shown in Figure 9. As 
in the two-channel ease study, the power spectral density peak magnitude 
increases for configurations with the higher (worse) rating, suggesting a 
proportional increase in pilot workload. 

It is noted that the peak power tends to occur at the dutch roll fre- 
quency for each configuration, indicating that this mode is clearly present 
if not dominant in the pilot's output. If this is the case this mode may 
be a contributing cause to the lateral PIO occurring for Configuration 3 
(refer to Table 1 for comments). 

The plots depicting noise contributions into the aileron and rudder 
deflection signals are shown in Figures 10 and 11. In addition to the 
large increase in peak spectral density of Configuration 3 over the other 
configurations, note that command disturbance noise is not dominant in the 
frequency range of maximum power as in the two-channel case (Figure 8). In 
the aileron deflection channel, pilot injected noise contribution exceeds 
the command disturbance noise contribution. This same trend is even more 
noticeable in the noise contribution plots for the rudder channel in Figure 
11, where the primary noise source is clearly pilot injected noise into the 
rudder channel . 

To summarize the data analysis of the identified models, there is evi- 
dence that the cause of the PIO and resultant poor pilot rating is self- 
induced coupling caused by rudder excitation of a dutch roll mode with 
level 3 flying qualities. Recall in the two-channel case study for Confi- 
guration 3 that no lateral PIO occurred when the rudder input was denied 
the pilot. The command disturbance in each case was identically provided 
using Equation (69). 

The frequency response of the pilot model, obtained from the approxi- 
mation of Equation (24), is shown for each configuration for the three- 
channel cases in Figures 12 and 13. Note that for for poorly rated Confi- 
guration 3 that pilot aileron deflection is out of phase at low frequencies 
with displayed bank angle error. 

As seen from the JIR analysis, the amount of information from the 
identification, validation, and analysis of models dbtained from actual 
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data sets is very large. Thus selectivity in analysis is essential, and 
the purpose of the modeling effort is paramount in this selection process. 

6^. CONCLUSIONS 

The fundamental conclusion from this research effort is that time 
series models and the analytical analysis tools they provide have the abil- 
ity to quantitatively evaluate pilot-in-the-loop situations by displaying 
key relationships affecting the stability and response of a multi-channel 
"piloted*' dynamic system. The NPD algorithm, in conjunction with the PAC 
selection rule, results in a parsimonious and stable multi-channel time 
series JIR model. This representation is unique if the existence condi- 
tions of Theorems 1 and 2 are met. Experimentally this requires sufficient 
and random-appearing excitation, physical delays in each path, and data 
realization sets which are stable. 

Analysis of case studies illustrated the application of the modeling 
process, and demonstrated how the dominant source of a lateral PIO may be 
identified using analysis tools presented in this paper. It is important 
to remember that the case study results were primarily intended to illus- 
trate the "application" of the identification process as opposed to a 
comprehensive evaluation of particular vehicle configurations. 

It is recommended that the joint innovations identification process be 
applied to a more varied data base, including actual flight test data and 
flight control system variations. Multi-channel applications which study 
manual control response of operators in training status may also be accom- 
plished. 
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8. APPENDIX 


PROOF OF THEOREM 1. Fro. Equation (10) „e have 


X(t) - 


I - G(M, z) 


-1 


fl(z)p(t) 


(A. 1 ) 


First the unique mapping between Equations (A.l) and (13) and (18) 

will be given, then the uniqueness conditions for the identified r(z) will 

be derived. Referring to Fi 8ure 1 and temporarily eliminating notation for 
arguments let 


K 1 ' (I ‘ G p G a G f> 
K 2 - « - G f G p G a ) 

K 3 = (I - G a G f G p ) 


(A. 2) 
(A. 3) 
(A. 4) 


Expand the subsystem blocks in Equation (18) to obtain 


X(t) - 


1*1, <*> 


11 

'21 


r„,( z ) 


r 3 1 ( z ) 


r i2'(z) r,,(z) 

r 22 (z) 

r ~->(z) r 0 ( z ) 


32 


13 

f 23 

r 33 


Pn^(z) 


p(t) 


(A. 5) 


Use direct substitution from Equation (A.l) and match entries with 
Equation (A. 5) to obtain 


F 11 = fl lfl H 

r !2 = VW 22 

r i3 = K 1 G p a 33 
r 21 " V®f0n 

r 22 \h a 22 

r 23 = K 2 G f G p fi 33 
r 31 “ V G . G f«u 
r 32 " V G « 22 

r 33 “ V *33 


(A. 6) 
(A. 7) 
(A. 8) 
(A. 9) 
(A. 10) 
(A. 11) 
(A. 12) 
(A. 13) 
(A. 14) 


nee 0 U are non-singular prewhitening filters, K singular implies 

" <B) 8in8Ular fr ° m EqUatl ° n (22) ‘ The reverse mapping is provided by the 
recursive relations in Equations (16) through (17). Note that if only two 

subsystems are present that Equations (A.6) through (A.14) yield the same 
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relationships given by Anderson‘S and Goto^. 

The final step in the proof is to show the uniqueness of T(z) and this 
will be done using the following result from Popov as communicated by 

Anderson^; 

For a nonsingular 

*«(«) = ^(z) (A. 15) 

there exists D(z) such that 

D*( z )D(z) - ^(z) (A. 16) 


det 


and there exists 


> 0 


(A. 17) 


^(w) = r(z) a p(0)r (z) 


= 


(A. 18) 


with T(z) and P(0) unique 


z = e 


r(z=») a i and P(0) > 0 


(A. 19) 


To apply this result to the JAR use the condition that there is a delay in 
every path , thus 


'ij 


z—°° 


= 0 


(A. 20) 


C. . =0 

ii z=« 


(A. 21) 


Substituting Equation (A. 20) into Equations (A. 6) through (A. 14), and sub- 
stituting Equation (A. 21) into Equations (11) and (12), we obtain 


ij 


Z=oo 


= 0 


(A. 22) 


ii 


z— 00 ^ 


(A. 23) 
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I Q ii j z=« " 1 (A. 24) 

Thus Equation (A. 19) is satisfied for the JAR. By the i.i.d. properties 

of p(t), P(0) is positive definite; and by the properties of a Toeplitz 
17 

Autocorrelation matrix 


satisfying Equation (A. 15). Therefore Popov's result applies and r(z) and 
P(0) are unique. Note that Anderson 10 has also shown that the block diago- 
nals of r (z) must be nonsingular. 

PROOF OF THEOREM 2. To prove that the joint process must be full rank 
for unique identification use Equation (28) 


xx 


r(z) a p(0)r (z) 


together with 


z=e 


jwA 


R xx (0) = E x <t)X T (t) 


r(z) A p(0)r*(z) 


If X(t) is less than full rank then a singular P(0) is implied. A singular 

P(0) makes one or more blocks of T(z) arbitrary. 

To prove the unique factorization is a necessary condition for unique 

identification use Equation (22) and the fact that P(0) is positive defin- 

23 

ite. Then there is a unitary transformation such that for some diagonal 
A(0) 


P(0) = UA(0)U T , A(0) > 0 


Therefore 


$ (bi) = 

xx 


r(z)UAA(0)u T r*(z)| 




1 z=e 


If P(0) is not diagonal, then the identified T(z) is 


f(z) - T(z) U 


(A. 25) 


(A. 26) 


(A. 27) 


where unitary matrix "U" depends on the correlation in P(0), and thus may 
not be unique. If P(0), however, is diagonal then P(0). = A(0), U = I, and 


21 



Hm r = r(z) u = r(z) 

M,N-h» 


(A. 28) 


Thus if the unitary matrix "U" is identity then a sufficient condition 

exists for the factorization to be unique. The "physically realizable" 

normalized minimum phase stable factor results as defined by Anderson 10 . 

_9. REFERENCES 

1. Onstott E., Faulkner, "Prediction, Evaluation, and Specification of 
Closed Loop and Multiaxis Flying Qualities" ,AFFDL-TR-78-3 , Wright- 
Patterson AFB.Ohio, Feb 1978. 

2. Frosch, R. /'Robots and People", Aeronautics and Astronautics, Vol 21, 
No. 7, pp. 34-36, Aug. 1983. 

3. Rouse W., Systems Engineering Models of Human-Machine Interaction, 
North Holland, 1980. 

4. Shinners S., "Modeling of Human Operator Performance Utilizing Time 
Series Analysis ".IEEE Tansactions on SMC, Vol SMC-4,No 5, pp 446- 
458, Sept 74. 

5. Goto N. ,"A Statistical Method Applied to Pilot Behavior Analysis in 
Multiloop Systems, "AIAA Journal of Guidance and Control, Vol 3, No 1, 
Jan-Feb 80, pp 62-68. 

6. Willsky A. ."Relationships Between Digital Signal Processing and Con- 
trol and Estimation Theory ".Proceedings of the IEEE, Vol 66, No 9, Sept 
78. 

7. Agarwal G., Gottlieb G., Osafo-Charles F., O'Niell W., "Application of 
Time-Series Modeling to Human Operator Dynamics", IEEE Transactions on 
Systems .Man, and Cybernetics, Vol SMC-10,No 12,pp849-860, Dec 1980. 

8. Akaike, H.,"0n the Use of a Linear Model for the Identification of 
Feedback Systems", Ann. Inst. Statist. Math. ,, Vol. 20, 1968, pp 425-439. 

9. Whittle P. , On the Fitting of Multivariate Regressions, and the Approx- 
imate Canonical Factorization of the Spectral Density Matrix", 
Blometrika.Vol 50,ppl29-134,1963. 

10. Anderson B.D.O., Ng T.S., Goodwin G.C., "Identif lability of MIMO 
Linear Dynamic Systems Operating in the Closed Loop", Automatica, Vol 
13, 1977, PP 477-485. 

11. Wiggins R. .Robinson E. /’Recursive Solution to the Multichannel Filter- 
ing Problem", Journal of Geophysical Research, Vol. 70, 1965, pp. 


22 



1885-1891. 

12. Morf M. , Vieira A., Lee D., Kailath T. ."Recursive Multi Channel Maximum 
Entropy Spectral Estimation" , IEEE Transactions on Geo Electronics, Vol 
GE-16,pp 85-94, Apr 78. 

13. Gelb A., Applied Optimal Estimation, MIT Press .Cambridge, Mass. ,1974 . 

14. Caines P.,Chan C., "Feedback Between Stationary Stochastic Processes", 
IEEE Transactions on Automatic Control, Vol AC-20, Aug 75, pp 498-508. 

15. Vorchik V., Fetisov V. .Shteinberg E ., "Identification of a Closed-loop 
Stochastic System",Automatika,No 7,pp41-52,July73. 

16. Anderson N.,"0n the Calculation of Filter Coefficients for Maximum 

Entropy Spectral Analysis", Geophysics, Vol. 39, No. 1, Feb 74, pp 
69-?72 . 

17. Kailath T.,"A View of Three Decades of Linear Filtering Theory", IEEE 
Transactions on Information Theory, Vol IT-20, ppl45-181 , Mar 1974. 

18. Morf M. , Vieira A., Kailath T. ."Covariance Characterization by Partial 
Autocorrelation Matrices", Annals of Statistics, Vol 6, No 3, pp 643- 
648, May 1978. 

19. Akaike, H. , "Canonical Correlation Analysis of Time Series and the Use 
of an Information Criterion", System Identification: Advances and Case 
Studies, Mehra ed.. Academic Press, 1976. 

20. Kashyap, "Inconsistency of Akaike Rule", IEEE Transactions on Automatic 
Control, Vol AC-25, No 5, Oct 1980, pp 997. 

21. "Military Specification-Flying Qualities of Piloted Airplanes", MIL- 
F-8785C, 1980. 

22. Schmidt D.K., "On the Use of the OCM Objective Function as a Pilot 
Rating Metric", 17th Annual Conference on Manual Control, UCLA, Los 
Angeles, June, 1981. 

23. Strang G. , Linear Algebra, Academic Press, 1978. 


23 



Table I. Three-channel case study configurations 


Configuration 
(Level ) 

Dutch Roll 
Parameters 

PR 

— — \ 

Comments 

£ 

cu N 

1: two-ch 

0.4 

2.02 

2 

"Responsive and predictable" 

1: three-ch 

— 

— 

4 

"Some coupling from rudder in to aileron axis, 
but mostly well behaved" 

2: two-ch 

0.1 

2.02 

HR 

"Some oscillations and overshoots when 
aggressive" 

2: three-ch 

“ “ 

“ 

7 

"Coupled overshoots between rudder and aileron 
and bank angle when aggressive; unpredictable 
and oscillatory bank angle made worse when 
aggressive on rudder." 

3: two-ch 

.02 

4.0* 

6 

"Overshoots and residual oscillations;" 
"unpredictable;" complex aileron inputs 
required for control" 

3: three-ch 

— 

— 

9 

"Closed loop unstable for task;" "excessive 
lateral PIO." 


*This is the frequency of the lateral PIO 






Table Z Two-channel case study results summary 


t'O 

Ln 


Configuration 
(Level ) 


PR 

Bank Angle 
Error 
. Variance 

deg 2 

Ail. Def. 

Variance 

, 2 
deg 

Cross 

Covariance 
Bank Angle 
to Aileron 
Deflection 

deg 2 

Maximum 
PSD Value 
of Aileron 
Defl ection 
2 

deg /rad/sec 

2 

20.9 

5.3 

8.2 

2.2 

4-5 

32.1 

6.6 

11.4 

2.9 

6 

77.1 

11.4 

20.8 

3.8 











Table 3 Covariance matrix summary and comparison 


Covariance * 

Configuration Matrix Order 

(Level) R vy (0) M 


5.3, 8.2 


8.2, 20.9 


" 6 . 6 , 11.4 
11.4, 32.1 


11.4, 20.8' 
20.8, 77.1 



Th ree-channel 



Covariance Matrix R (0) 

xx v ' 

[equation (4.43)] 


4.93, 8.55, 8.83 
8.55, 53.44, 23.96 
8.83, 23.96, 25.0 


'19.71, 24.3 , 22.46' 
24.3 , 138.1 , 43.87 
22.46, 43.87, 34.17 


11.69, 15.62, 8.21 
15.62, 541.1 , 41.78 
8.21, 41.78, 30.5 


Normalized P(0) Matrix 
[equation (4.44)] 


0.10, 0.15 


0.08, 1.0 J 


0.11, 0.15' 
1.0 , 0.11 
0.11, 1.0 


- .19 , .0 

1.0 , .0 
.011, 1.0 


order obtained for PAC selection rule 






Figure 1 Multi-channel piloted closed-loop system model 

Figure 2 Multi-channel lateral axis tracking display 

Figure 3 Order selection rule 

Figure 4 Rudder channel actual vs. model output: 3-ch case 

study configuration 3 

Figure 5 Autocorrelation matrix vs. lag: 2-ch case study 

configuration 3 

Figure 6 Residual autocorrelation matrix vs. lag: 2-ch case 

study configuration 3 

Figure 7 Total aileron deflection power: 2-ch case study 

Figure 8 Noise contribution to aileron deflection PSC: 2-ch 

case study configuration 3 

Figure 9 Total aileron deflection power: 3-ch case study 

Figure IQ Noise contribution to aileron deflection PSD: 3-ch 

case study configuration 3 

Figure 11 Fjoise contribution to rudder deflection PSD: 3-ch 

case study configuration 3 

Figure 12 Frequency response magnitude 6 a /e g : 3-ch case study 

Figure 13 Frequency response phase 6 /e 3 : 3-ch case study 

a a 
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